7-4 Days alive and free of hospital by day 28

This outcome was defined as the number of days alive and free of hospital from randomisation to 28 days, calculated as 28 minus the number of days spent in hospital. All patients dying within 28 days will be assigned zero free days.

Authors
Affiliations
James Totterdell

University of Sydney

Rob Mahar

University of Melbourne

Published

July 14, 2022

Preamble

Load packages
library(tidyverse)
library(labelled)
library(kableExtra)
library(cmdstanr)
library(posterior)
library(bayestestR)
library(bayesplot)
library(matrixStats)
library(plotly)
library(lubridate)
library(ggdist)
library(patchwork)

theme_set(theme_classic(base_size = 10, base_family = "Palatino") +
  theme(panel.grid = element_blank(),
        strip.background = element_blank()))
bayesplot_theme_set(theme_set(theme_classic(base_size = 10, base_family = "Palatino") +
  theme(panel.grid = element_blank(),
        strip.background = element_blank())))

color_scheme_set("red")
options(digits = 4)
Prepare dataset
devtools::load_all()
all_dat <- read_all_no_daily()

acs_itt_dat <- all_dat %>% 
  filter_acs_itt() %>%
  transmute_model_cols_grp_aus_nz()
acs_itt_nona_dat <- acs_itt_dat %>% 
  filter(!is.na(out_dafh))

# Concurrent enrolments for C2
acs_itt_concurc2_dat <- acs_itt_dat %>%
  filter_concurrent_intermediate()
acs_itt_concurc2_nona_dat <- acs_itt_concurc2_dat %>% 
  filter(!is.na(out_dafh))

# Concurrent enrolments for C3
acs_itt_concurc3_dat <- acs_itt_dat %>%
  filter_concurrent_std_aspirin()
acs_itt_concurc3_nona_dat <- acs_itt_concurc3_dat %>% 
  filter(!is.na(out_dafh))

# Concurrent enrolments for C4
acs_itt_concurc4_dat <- acs_itt_dat %>%
  filter_concurrent_therapeutic()
acs_itt_concurc4_nona_dat <- acs_itt_concurc4_dat %>% 
  filter(!is.na(out_dafh))
Load models
ordmod0 <- cmdstan_model("stan/ordinal/logistic_cumulative.stan") # No epoch or site
ordmod_epoch <- cmdstan_model("stan/ordinal/logistic_cumulative_epoch.stan") # Epoch only
ordmod <- cmdstan_model("stan/ordinal/logistic_cumulative_epoch_site.stan") # Full model
Functions
make_summary_table <- function(dat, format = "html") {
  tdat <- dat %>%
  group_by(CAssignment = factor(CAssignment, 
                                levels = c("C1", "C2", "C3", "C4"),
                                labels = intervention_labels2()$CAssignment[-1])) %>%
  summarise(
    Patients = n(),
    Known = sum(!is.na(out_dafh)),
    Deaths = sprintf(
      "%i (%.0f%%)", sum(D28_death, na.rm = TRUE), 100 * mean(D28_death, na.rm = TRUE)),
    `DAFH, Median (Q1, Q3)` = sprintf(
      "%.0f (%.0f, %.0f)", 
      median(out_dafh, na.rm = T), 
      quantile(out_dafh, 0.25, na.rm = TRUE), 
      quantile(out_dafh, 0.75, na.rm = TRUE))
  ) %>%
  bind_rows(
    dat %>%
  group_by(CAssignment = "Overall") %>%
  summarise(
    Patients = n(),
    Known = sum(!is.na(out_dafh)),
    Deaths = sprintf(
      "%i (%.0f%%)", sum(D28_death, na.rm = TRUE), 100 * mean(D28_death, na.rm = TRUE)),
    `DAFH, Median (Q1, Q3)` = sprintf(
      "%.0f (%.0f, %.0f)", 
      median(out_dafh, na.rm = T), 
      quantile(out_dafh, 0.25, na.rm = TRUE), 
      quantile(out_dafh, 0.75, na.rm = TRUE))
  )
  ) %>%
  rename(`Anticoagulation\nintervention` = CAssignment)
  kable(
    tdat,
    format = format,
    align = "lrrrr",
    booktabs = TRUE,
    linesep = ""
  ) %>%
    kable_styling(
      font_size = 9,
      latex_options = "HOLD_position"
    ) %>%
    row_spec(nrow(tdat), bold = T)
}


make_primary_model_data <- function(
    dat, 
    vars = NULL,
    beta_sd_var = NULL,
    ctr = contr.orthonorm,
    p_mult = 2) {
  
  X <- make_X_design(dat, vars, ctr)
  attX <- attr(X, "contrasts")
  X <- X[, -1]
  attr(X, "contrasts") <- attX
  nXtrt <- sum(grepl("rand", colnames(X)))
  
  beta_sd <- c(rep(1, nXtrt), beta_sd_var)
  
  epoch <- dat$epoch
  M_epoch <- max(dat$epoch)
  region <- dat$ctry_num
  M_region <- max(region)
  site <- dat$site_num
  M_site <- max(site)
  region_by_site <- dat %>% 
    dplyr::count(ctry_num, site_num) %>% 
    pull(ctry_num)

  y_raw <- ordered(dat$out_dafh)
  y_mod <- as.integer(y_raw)
  N <- dim(X)[1]
  K <- dim(X)[2]
  unique_y <- length(unique(y_mod))
  list(N = N, K = K, J = unique_y, X = X, y = y_mod, y_raw = y_raw,
       M_region = M_region,
       M_site = M_site, site = site,
       M_epoch = M_epoch, epoch = epoch,
       region_by_site = region_by_site,
       p_par = p_mult * rep(1 / unique_y, unique_y),
       beta_sd = beta_sd)  
}


fit_primary_model <- function(
    dat, 
    ctr = contr.orthonorm,
    save = FALSE
  ) {
  
  mdat <- make_primary_model_data(
    dat, c("inelgc3", "agegte60", "ctry"), c(10, 2.5, 1, 1), ctr)
  snk <- capture.output(
    mfit <- ordmod$sample(
      data = mdat,
      seed = 813508,
      chains = 8,
      parallel_chains = min(8, parallel::detectCores()),
      iter_warmup = 1000,
      iter_sampling = 2500,
      refresh = 0,
      adapt_delta = 0.98
  ))
  if (save) {
    mfit$save_object(file.path("outputs", "models", "secondary", "7-4.rds"))  
  }
  mdrws <- as_draws_rvars(mfit$draws())
  
  # Add names
  epoch_map <- dat %>% dplyr::count(epoch, epoch_lab)
  site_map <- dat %>% dplyr::count(site_num, site)
  names(mdrws$beta) <- colnames(mdat$X)
  names(mdrws$gamma_epoch) <- epoch_map$epoch_lab
  names(mdrws$gamma_site) <- site_map$site
  
  # Convert to treatment log odds ratios
  mdrws$Acon <- attr(mdat$X, "contrasts")$randA %**% mdrws$beta[grepl("randA[0-9]", names(mdrws$beta))]
  mdrws$Ccon <- attr(mdat$X, "contrasts")$randC %**% mdrws$beta[grepl("randC[0-9]", names(mdrws$beta))]
  mdrws$trtA <- mdrws$Acon[-1] - mdrws$Acon[1]
  mdrws$trtC <- mdrws$Ccon[-1] - mdrws$Ccon[1]
  
  mdrws$compare <- c(
    "Intermediate vs low" = exp(mdrws$trtC[1]),
    "Low with aspirin vs low" = exp(mdrws$trtC[2]),
    "Therapeutic vs low" = exp(mdrws$trtC[3]),
    "Intermediate vs low with aspirin" = exp(mdrws$trtC[1] - mdrws$trtC[2]),
    "Intermediate vs therapeutic" = exp(mdrws$trtC[1] - mdrws$trtC[3]),
    "Low with aspirin vs therapeutic" = exp(mdrws$trtC[2] - mdrws$trtC[3])
  )
  mdrws$OR <- c(
    setNames(exp(mdrws$trtC), c("Intermediate", "Low with aspirin", "Therapeutic")),
    "Ineligible aspirin" = exp(mdrws$beta[which(grepl("inelg", colnames(mdat$X)))]),
    "Age \u2265 60" = exp(mdrws$beta[which(grepl("age", colnames(mdat$X)))]),
    setNames(exp(mdrws$beta[which(grepl("ctry", colnames(mdat$X)))]), c("AU/NZ", "NP"))
  )
  return(list(dat = mdat, fit = mfit, drws = mdrws))
}

Outcome Derivation

The outcome is calculated for a patient as:

  • missing, if day 28 mortality is unknown (D28_PatientStatus) or if the patient was known to have been alive at day 28 but number of days in hospital (D28_OutcomeTotalDaysHospitalised) is unknown
  • 0 if they died by day 28 (D28_PatientStatus)
  • 28 - min(28, D28_OutcomeTotalDaysHospitalised) otherwise

Descriptive Analyses

Overall

As a cross check, Figure 1 plots the number of days hospitalised (as reported in D28_OutcomeTotalDaysHospitalised) against the day of discharge from the index admission.

Compare D28_OutcomeTotalDaysHospitalised to days hospitalised per dischage date.
pdat <- acs_itt_dat %>%
  filter(D28_death != 1) %>%
  dplyr::count(D28_day = pmin(28, D28_OutcomeTotalDaysHospitalised), DIS_day = pmin(28, DIS_day))
ggplot(pdat, aes(D28_day, DIS_day)) +
  geom_point(aes(size = n), shape = 21) +
  geom_abline() +
  scale_x_continuous("Number of days hospitalised (D28_OutcomeTotalDaysHospitalised)",
                     breaks = seq(0, 30, 2)) +
  scale_y_continuous("Day of discharge from index admission",
                     breaks = seq(0, 30, 2)) +
  labs(size = "Count")

The overall distribution of days alive and free of hospital (DAFH) are reported in Figure 2 for all participants in the ACS-ITT set.

Overall distribution of DAFH
pdat <- acs_itt_nona_dat %>%
  dplyr::count(dafh = ordered(out_dafh, levels = 0:27), .drop = F) %>%
  mutate(p = n / sum(n))
p <- ggplot(pdat, aes(dafh, p)) +
  geom_col() +
  labs(
    x = "Days alive and free of hospital", 
    y = "Proportion"
  )
path <- file.path("outputs", "figures", "outcomes", "secondary")
ggsave(file.path(path, "outcome-7-4-dist-overall.pdf"), p, height = 2.5, width = 6)
p

Summary of DAFH outcome by arm
save_tex_table(
  make_summary_table(acs_itt_dat, "latex"), 
  file.path("outcomes", "secondary", "7-4-anticoagulation-summary"))
make_summary_table(acs_itt_dat)
Anticoagulation intervention Patients Known Deaths DAFH, Median (Q1, Q3)
Low-dose 610 596 19 (3%) 23 (21, 24)
Intermediate-dose 613 603 15 (2%) 23 (21, 24)
Low-dose with aspirin 283 280 10 (4%) 22 (20, 24)
Therapeutic-dose 50 50 6 (12%) 22 (19, 24)
Overall 1556 1529 50 (3%) 23 (21, 24)
Table 1: Summary of days alive and free of hospital to day 28 by treatment group.

Age

DAFH by age
pdat <- acs_itt_nona_dat %>%
  dplyr::count(
    agegrp = cut(AgeAtEntry, c(18, seq(25, 75, 5), 100), include.lowest = T, right = F),
    dafh = fct_rev(ordered(out_dafh, levels = 0:27)), 
    .drop = F) %>%
  group_by(agegrp) %>%
  mutate(p = n / sum(n))
pdat2 <- pdat %>%
  group_by(agegrp) %>%
  summarise(n = sum(n))
p1 <- ggplot(pdat2, aes(agegrp, n)) +
  geom_col(colour = "grey40", fill = "grey40") +
  geom_vline(xintercept = 60, linetype = 2) +
  labs(y = "Number of\nparticipants",
       x = "Age at entry") +
  geom_vline(xintercept = 8.5, linetype = 2) +
  theme(axis.text.x = element_text(hjust = 1, angle = 45))
p2 <- ggplot(pdat, aes(agegrp, p)) +
  geom_col(aes(fill = dafh)) +
  labs(x = "Age", y = "Cumulative\nproportion") +
  scale_fill_viridis_d("DAFH", option = "A", direction = -1) +
  guides(fill = guide_legend(reverse = TRUE)) +
  geom_vline(xintercept = 8.5, linetype = 2) +
  theme(axis.text.x = element_text(hjust = 1, angle = 45)) +
  theme(legend.key.size = unit(0.25, "lines"))
p <- p1 | p2
path <- file.path("outputs", "figures", "outcomes", "secondary")
ggsave(file.path(path, "7-4-age.pdf"), p, height = 2.5, width = 6)
p

Country

DAFH by country
pdat <- acs_itt_nona_dat %>%
  dplyr::count(
    country,
    dafh = ordered(out_dafh, levels = 0:27), 
    .drop = F) %>%
  group_by(country) %>%
  mutate(p = n / sum(n)) %>%
  mutate(cp = cumsum(p)) %>%
  ungroup()
pdat2 <- pdat %>%
  group_by(country) %>%
  summarise(n = sum(n))

p1 <- ggplot(pdat2, aes(country, n)) +
  geom_col() +
    labs(
      y = "Number of\nparticipants", 
      x = "Country of enrolment")

p2 <- ggplot(pdat, aes(country, p)) +
  geom_col(aes(fill = fct_rev(dafh))) +
  labs(x = "Anticoagulation intervention", y = "Cumulative\nproportion") +
  scale_fill_viridis_d("DAFH", option = "A", direction = -1) +
  guides(fill = guide_legend(reverse = TRUE)) +
  theme(legend.key.size = unit(0.25, "lines"))
p <- p1 | p2
path <- file.path("outputs", "figures", "outcomes", "secondary")
ggsave(file.path(path, "7-4-country.pdf"), p, height = 2.5, width = 6)
p

Site

DAFH by site
pdat <- all_dat %>%
  filter_acs_itt() %>%
  filter(!is.na(out_dafh)) %>%
  dplyr::count(
    Country = factor(PT_CountryName, levels = c("India", "Australia", "Nepal", "New Zealand"),
                     labels = c("India", "Australia", "Nepal", "New\nZealand")),
    Site = fct_infreq(Location),
    dafh = ordered(out_dafh, levels = 0:27)) %>%
  complete(dafh = ordered(0:27), nesting(Country, Site), fill = list(n = 0)) %>%
  group_by(Country, Site) %>%
  mutate(p = n / sum(n)) %>%
  mutate(cp = cumsum(p)) %>%
  ungroup() %>%
  mutate(
    Country = droplevels(Country),
    Site = droplevels(Site)
  )
pdat2 <- pdat %>%
  group_by(Country, Site) %>%
  summarise(n = sum(n)) %>%
  ungroup()
p1 <- ggplot(pdat2, aes(Site, n)) +
  facet_grid( ~ Country, scales = "free_x", space = "free_x") +
  geom_col() +
    labs(
      y = "Number of\nparticipants", 
      x = "") +
  theme(axis.text.x = element_text(angle = 45, hjust = 1),
        panel.border = element_rect(fill = NA))
p2 <- ggplot(pdat, aes(Site, p)) +
  facet_grid( ~ Country, scales = "free_x", space = "free_x") +
  geom_col(aes(fill = fct_rev(dafh))) +
  labs(x = "Anticoagulation intervention", y = "Cumulative\nproportion") +
  scale_fill_viridis_d("DAFH", option = "A", direction = -1) +
  guides(fill = guide_legend(reverse = TRUE, ncol = 1)) +
  theme(axis.text.x = element_text(hjust = 1, angle = 45)) +
  theme(legend.key.size = unit(0.25, "lines"))
p <- p1 / p2 +
  plot_layout(guides = 'collect')
path <- file.path("outputs", "figures", "outcomes", "secondary")
ggsave(file.path(path, "7-4-country-site.pdf"), p, height = 4, width = 6.25)
p

Calendar Time

DAFH by calendar date
pdat <- acs_itt_nona_dat %>%
  dplyr::count(
    yr = year(RandDate), mth = month(RandDate),
    dafh = ordered(out_dafh, levels = 0:27), 
    .drop = F) %>%
  group_by(yr, mth) %>%
  mutate(p = n / sum(n)) %>%
  mutate(cp = cumsum(p)) %>%
  ungroup()
p1 <- pdat %>%
  group_by(yr, mth) %>%
  summarise(n = sum(n)) %>%
  ggplot(., aes(mth, n))  +
  facet_grid( ~ yr, drop = T, scales = "free_x", space = "free") +
    geom_col() +
    labs(
      y = "Number of\nparticipants", 
      x = "Calendar date (month of year)") +
  scale_x_continuous(breaks = 1:12)
p2 <- ggplot(pdat, aes(mth, p)) +
  facet_grid( ~ yr, drop = T, scales = "free_x", space = "free") +
  geom_col(aes(fill = fct_rev(dafh))) +
  labs(x = "Calendar date (month of year)", y = "Cumulative\nproportion") +
  scale_fill_viridis_d("DAFH", option = "A", direction = -1) +
  guides(fill = guide_legend(reverse = TRUE)) +
  theme(legend.key.size = unit(0.25, "lines")) +
  scale_x_continuous(breaks = 1:12)
p <- p1 | p2
path <- file.path("outputs", "figures", "outcomes", "secondary")
ggsave(file.path(path, "7-4-calendar-time.pdf"), p, height = 2, width = 6)
p

Anti-coagulation

DAFH by intervention
pdat <- acs_itt_nona_dat %>%
  dplyr::count(
    CAssignment = factor(CAssignment, labels = str_replace(intervention_labels()$CAssignment[-1], "<br>", "\n")),
    dafh = ordered(out_dafh, levels = 0:27), 
    .drop = F) %>%
  group_by(CAssignment) %>%
  mutate(p = n / sum(n)) %>%
  mutate(cp = cumsum(p)) %>%
  ungroup()
p <- ggplot(pdat, aes(CAssignment, p)) +
  geom_col(aes(fill = fct_rev(dafh))) +
  labs(x = "Anticoagulation intervention", y = "Cumulative proportion") +
  scale_fill_viridis_d("Days alive and\nfree of hospital", option = "A", direction = -1) +
  guides(fill = guide_legend(reverse = TRUE)) +
  theme(legend.key.size = unit(0.75, "lines"))
pth <- file.path("outputs", "figures", "outcomes", "secondary")
ggsave(file.path(pth, "outcome-7-4-descriptive.pdf"), p, height = 3, width = 6)
p

DAFH by intervention
p <- ggplot(pdat, aes(dafh, cp)) +
  geom_step(aes(colour = CAssignment, group = CAssignment)) +
  labs(x = "Days alive and free of hospital", y = "Cumulative proportion") +
  scale_colour_viridis_d(
    "Days alive and\nfree of hospital", option = "A", begin = 0, end = 0.8) +
  guides(fill = guide_legend(reverse = TRUE))
pth <- file.path("outputs", "figures", "outcomes", "secondary")
ggsave(file.path(pth, "outcome-7-4-descriptive2.pdf"), p, height = 3, width = 6)
p

Sample Cumulative Logits

Some evidence of a shift from proportional odds at higher values of DAFH.

Plot sample cumulative logits
trt_counts <- acs_itt_nona_dat %>%
  dplyr::count(CAssignment, out_dafh) %>%
  complete(CAssignment, out_dafh, fill = list(n = 0)) %>%
  group_by(CAssignment) %>%
  mutate(p = n / sum(n))
trt_logit <- trt_counts %>% 
  group_by(CAssignment) %>% 
  mutate(clogit = logit(cumsum(p))) %>%
  group_by(out_dafh) %>%
  mutate(rel_clogit = clogit - mean(clogit)) %>%
  filter(out_dafh != 27)

ggplot(trt_logit, aes(CAssignment, rel_clogit)) +
  facet_wrap( ~ out_dafh) +
  geom_point() +
  labs(y = "Relative (to mean) sample cumulative logit")

Plot sample cumulative logits
ggplot(trt_logit, aes(out_dafh, rel_clogit)) +
  facet_wrap( ~ CAssignment) +
  geom_point() +
  labs(y = "Relative (to mean) sample cumulative logit")

Modelling

The primary analysis uses the ACS-ITT set, so all participants have been randomised to the anticoagulation domain and some participants may have been randomised to the antiviral domain.

Primary Model

Fit primary model
res <- fit_primary_model(acs_itt_nona_dat, save = TRUE)
Odds ratio summary table
save_tex_table(
  odds_ratio_summary_table_rev(res$drws$OR, "latex"),
  "outcomes/secondary/7-4-primary-model-acs-itt-summary-table")
odds_ratio_summary_table_rev(res$drws$OR)
Parameter Median 95% CrI Mean (SD) Pr(OR > 1)
Intermediate 1.14 (0.94, 1.39) 1.15 (0.12) 0.90
Low with aspirin 1.05 (0.80, 1.39) 1.06 (0.15) 0.65
Therapeutic 0.74 (0.41, 1.35) 0.77 (0.24) 0.16
Ineligible aspirin 1.08 (0.56, 2.05) 1.14 (0.39) 0.59
Age ≥ 60 0.58 (0.47, 0.71) 0.59 (0.06) 0.00
AU/NZ 0.75 (0.33, 1.73) 0.82 (0.37) 0.25
NP 0.68 (0.23, 2.62) 0.86 (0.72) 0.26
Posterior summaries for model parameters (fixed-effects).
Odds ratio summary for epoch and site
p <- plot_epoch_site_terms(
  exp(res$drws$gamma_epoch),
  exp(res$drws$gamma_site),
  factor(res$dat$region_by_site, 
         labels = c("India", "Australia\nNew Zealand", "Nepal"))
)
pth <- file.path("outputs", "figures", "outcomes", "secondary", "7-4-primary-model-epoch-site-terms.pdf")
ggsave(pth, p, width = 6, height = 4.5)
p

Posterior contrasts
plot_or_densities(res$drws$compare)

Code
res$fit$summary(
  c("alpha", "beta", "gamma_epoch", "gamma_site", "tau_epoch", "tau_site")) %>%
  print(n = Inf)
# A tibble: 76 × 10
   variable    mean   median    sd   mad      q5     q95  rhat ess_bulk ess_tail
   <chr>      <dbl>    <dbl> <dbl> <dbl>   <dbl>   <dbl> <dbl>    <dbl>    <dbl>
 1 alpha[1] -4.03   -4.03    0.306 0.301 -4.52   -3.52    1.00    2367.    4074.
 2 alpha[2] -4.01   -4.02    0.306 0.302 -4.50   -3.50    1.00    2355.    4075.
 3 alpha[3] -4.00   -4.00    0.306 0.302 -4.48   -3.49    1.00    2351.    4140.
 4 alpha[4] -3.98   -3.98    0.305 0.302 -4.47   -3.47    1.00    2347.    4113.
 5 alpha[5] -3.90   -3.90    0.304 0.299 -4.38   -3.39    1.00    2310.    4173.
 6 alpha[6] -3.88   -3.89    0.304 0.300 -4.37   -3.38    1.00    2302.    4129.
 7 alpha[7] -3.81   -3.81    0.302 0.297 -4.29   -3.31    1.00    2275.    4082.
 8 alpha[8] -3.75   -3.76    0.301 0.297 -4.24   -3.25    1.00    2260.    4008.
 9 alpha[9] -3.66   -3.67    0.300 0.294 -4.14   -3.16    1.00    2220.    3986.
10 alpha[1… -3.64   -3.64    0.299 0.293 -4.11   -3.14    1.00    2209.    3928.
11 alpha[1… -3.52   -3.53    0.297 0.293 -3.99   -3.03    1.00    2175.    3924.
12 alpha[1… -3.48   -3.48    0.296 0.292 -3.95   -2.98    1.00    2163.    3860.
13 alpha[1… -3.43   -3.44    0.296 0.290 -3.91   -2.94    1.00    2148.    3810.
14 alpha[1… -3.35   -3.36    0.294 0.290 -3.82   -2.86    1.00    2125.    3791.
15 alpha[1… -3.13   -3.14    0.292 0.287 -3.59   -2.64    1.00    2111.    3711.
16 alpha[1… -2.99   -2.99    0.290 0.285 -3.45   -2.50    1.00    2083.    3672.
17 alpha[1… -2.83   -2.83    0.289 0.283 -3.28   -2.34    1.00    2065.    3703.
18 alpha[1… -2.50   -2.50    0.286 0.278 -2.95   -2.02    1.00    2047.    3776.
19 alpha[1… -2.05   -2.06    0.284 0.276 -2.50   -1.57    1.00    2050.    3709.
20 alpha[2… -1.50   -1.51    0.282 0.277 -1.95   -1.03    1.00    2034.    3704.
21 alpha[2… -0.781  -0.787   0.280 0.274 -1.23   -0.310   1.00    2034.    3735.
22 alpha[2…  0.231   0.226   0.280 0.275 -0.213   0.701   1.00    2060.    3742.
23 alpha[2…  1.14    1.14    0.283 0.275  0.695   1.62    1.00    2085.    3745.
24 alpha[2…  2.76    2.75    0.300 0.295  2.28    3.26    1.00    2315.    4362.
25 alpha[2…  5.30    5.27    0.502 0.494  4.51    6.16    1.00    5648.    9034.
26 beta[1]  -0.259  -0.258   0.233 0.236 -0.639   0.124   1.00   14764.   15129.
27 beta[2]   0.0348  0.0351  0.111 0.110 -0.147   0.218   1.00   13839.   14787.
28 beta[3]   0.202   0.203   0.136 0.135 -0.0212  0.425   1.00   11940.   13893.
29 beta[4]  -0.174  -0.173   0.405 0.405 -0.838   0.485   1.00   15857.   14629.
30 beta[5]  -0.105  -0.104   0.243 0.244 -0.505   0.295   1.00   28996.   15863.
31 beta[6]   0.0750  0.0760  0.334 0.337 -0.476   0.617   1.00   29688.   16193.
32 beta[7]  -0.541  -0.541   0.105 0.105 -0.715  -0.371   1.00   22612.   14452.
33 beta[8]  -0.286  -0.294   0.423 0.423 -0.971   0.415   1.00    6844.   10122.
34 beta[9]  -0.350  -0.383   0.612 0.570 -1.29    0.693   1.00    9772.   12454.
35 gamma_e…  0       0       0     0      0       0      NA         NA       NA 
36 gamma_e… -0.218  -0.216   0.209 0.203 -0.565   0.122   1.00   13141.   14932.
37 gamma_e… -0.386  -0.384   0.216 0.212 -0.744  -0.0325  1.00   12306.   15291.
38 gamma_e… -0.202  -0.203   0.209 0.209 -0.543   0.143   1.00   12854.   15539.
39 gamma_e… -0.339  -0.339   0.196 0.196 -0.661  -0.0193  1.00   12483.   15628.
40 gamma_e… -0.537  -0.536   0.207 0.205 -0.884  -0.200   1.00   12284.   16486.
41 gamma_e… -0.614  -0.614   0.215 0.214 -0.968  -0.263   1.00   11700.   14932.
42 gamma_e… -0.699  -0.698   0.209 0.207 -1.04   -0.355   1.00   11636.   14836.
43 gamma_e… -1.38   -1.38    0.224 0.224 -1.75   -1.01    1.00   10313.   15333.
44 gamma_e… -1.56   -1.56    0.224 0.225 -1.93   -1.19    1.00   10098.   14636.
45 gamma_e… -0.969  -0.972   0.245 0.247 -1.37   -0.564   1.00   12096.   16271.
46 gamma_e… -0.728  -0.728   0.309 0.307 -1.24   -0.216   1.00   14522.   16103.
47 gamma_s… -0.434  -0.441   0.287 0.277 -0.894   0.0465  1.00    4269.    6968.
48 gamma_s…  0.0703  0.0620  0.272 0.262 -0.358   0.525   1.00    3606.    5699.
49 gamma_s… -0.310  -0.310   0.392 0.379 -0.958   0.327   1.00    9072.   11499.
50 gamma_s… -0.198  -0.205   0.271 0.258 -0.628   0.263   1.00    3452.    6273.
51 gamma_s…  0.192   0.181   0.243 0.231 -0.188   0.604   1.00    2602.    4698.
52 gamma_s… -0.818  -0.817   0.305 0.301 -1.32   -0.316   1.00    5343.    8908.
53 gamma_s…  1.19    1.18    0.274 0.263  0.765   1.66    1.00    3186.    5573.
54 gamma_s…  0.462   0.454   0.317 0.308 -0.0379  1.00    1.00    5139.    7986.
55 gamma_s…  0.723   0.710   0.273 0.262  0.297   1.19    1.00    3645.    5908.
56 gamma_s… -0.450  -0.454   0.324 0.315 -0.976   0.0874  1.00    5700.    8779.
57 gamma_s…  0.592   0.582   0.300 0.292  0.119   1.10    1.00    4694.    6890.
58 gamma_s… -0.252  -0.135   0.392 0.258 -1.03    0.192   1.00   14816.   16811.
59 gamma_s…  0.0525  0.0180  0.285 0.196 -0.385   0.568   1.00   24551.   17235.
60 gamma_s…  0.0174  0.00269 0.321 0.206 -0.497   0.574   1.00   26417.   16806.
61 gamma_s…  0.0408  0.0139  0.297 0.199 -0.421   0.562   1.00   27477.   17553.
62 gamma_s… -0.135  -0.0593  0.333 0.215 -0.770   0.308   1.00   20678.   16621.
63 gamma_s…  0.166   0.0908  0.309 0.217 -0.230   0.769   1.00   18778.   17195.
64 gamma_s… -0.120  -0.0539  0.324 0.214 -0.731   0.326   1.00   22194.   17162.
65 gamma_s…  0.0885  0.0325  0.347 0.211 -0.402   0.730   1.00   22614.   16515.
66 gamma_s… -0.187  -0.0929  0.349 0.231 -0.865   0.244   1.00   18447.   17081.
67 gamma_s…  0.0976  0.0376  0.322 0.206 -0.356   0.699   1.00   25117.   17831.
68 gamma_s…  0.150   0.0604  0.367 0.228 -0.324   0.868   1.00   20749.   16890.
69 gamma_s…  0.106   0.0613  0.249 0.188 -0.251   0.559   1.00   20773.   18295.
70 gamma_s… -0.0677 -0.0249  0.314 0.207 -0.639   0.409   1.00   23874.   17451.
71 gamma_s… -0.739  -0.683   0.623 0.564 -1.85    0.165   1.00   10772.   12214.
72 gamma_s…  0.405   0.407   0.615 0.533 -0.626   1.41    1.00   12780.   13080.
73 tau_epo…  0.408   0.390   0.124 0.110  0.245   0.635   1.00    7665.   11712.
74 tau_sit…  0.710   0.679   0.193 0.170  0.459   1.07    1.00    5315.    9463.
75 tau_sit…  0.316   0.276   0.229 0.229  0.0271  0.749   1.00    7860.    8643.
76 tau_sit…  1.03    0.884   0.649 0.476  0.338   2.23    1.00   10934.   10619.
Code
res$fit$diagnostic_summary()
$num_divergent
[1] 0 0 0 0 0 0 0 0

$num_max_treedepth
[1] 0 0 0 0 0 0 0 0

$ebfmi
[1] 0.8161 0.8034 0.7978 0.8037 0.7995 0.7954 0.8154 0.8189
Code
mcmc_trace(res$drws["beta"])

Code
mcmc_trace(res$drws["alpha"])

Code
mcmc_trace(res$drws["gamma_site"])

Code
mcmc_trace(res$drws["gamma_epoch"])

Posterior Predictive

Code
y_raw <- as.integer(levels(res$dat$y_raw))
y_ppc <- res$drws$y_ppc
y_ppc_raw <- rfun(\(x) y_raw[x])(y_ppc)
ppc_dat <- bind_cols(acs_itt_nona_dat, tibble(y_ppc = y_ppc_raw)) %>%
  mutate(CAssignment = factor(
    CAssignment, 
    levels = names(intervention_labels_short_break()$CAssignment),
    labels = intervention_labels_short_break()$CAssignment))

grp_ppc <- function(grp, d = 0:27) {
  ppc_dat %>%
    group_by(grp = {{grp}}) %>%
    summarise(
      y_lteq = map_dbl(d, ~ mean(out_dafh <= .x)),
      ypp_lteq = map(d, ~ rvar_mean(y_ppc <= .x)),
      y_eq = map_dbl(d, ~ mean(out_dafh == .x)),
      ypp_eq = map(d, ~ rvar_mean(y_ppc == .x))
    ) %>%
    unnest(c(ypp_lteq, ypp_eq)) %>%
    mutate(days = d) %>%
    ungroup() %>%
    pivot_longer(
      y_lteq:ypp_eq, 
      names_to = c("response", "event"),
      names_sep = "_",
      values_to = "posterior")
}

plot_grp_ppc <- function(dat, lab = "", xlab = "Probability") {
  ggplot(dat %>% 
           filter(response == "ypp", event == "eq"), 
         aes(y = days)) +
    facet_wrap( ~ grp, nrow = 1, scales = "free_x") +
    stat_slabinterval(aes(xdist = posterior), fatten_point = 1)  +
    geom_point(data = dat %>% filter(response == "y", event == "eq"), 
               aes(y = days, x = mean(posterior)),
               colour = "red",
               shape = 23) +
    scale_x_continuous(
      xlab, breaks = c(0, 0.3, 0.6),
      sec.axis = sec_axis(
        ~ . , name = lab, breaks = NULL, labels = NULL)) +
    scale_y_continuous("Days") +
    theme(strip.text = element_text(size = rel(0.7)),
          axis.title.x = element_text(size = rel(0.7)),
          axis.text.x = element_text(size = rel(0.65)),
          axis.title.y = element_text(size = rel(0.75)),
          axis.title.x.bottom = element_blank())
}

pp_C <- grp_ppc(CAssignment)
pp_ctry <- grp_ppc(country)
pp_epoch <- grp_ppc(epoch)
pp_site <- grp_ppc(site) %>%
  left_join(ppc_dat %>% dplyr::count(site, country), 
            by = c("grp" = "site"))

p1 <- plot_grp_ppc(pp_C, "Anticoagulation", "")
p2 <- plot_grp_ppc(pp_ctry, "Country", "") 
p3 <- plot_grp_ppc(pp_epoch, "Epoch", "")
p4 <- plot_grp_ppc(
  pp_site %>% filter(country == "IN"), "Sites India", "")
p5 <- plot_grp_ppc(
  pp_site %>% filter(country == "AU"), "Sites Australia", "")
p6 <- plot_grp_ppc(
  pp_site %>% filter(country == "NP"), "Sites Nepal", "")
p7 <- plot_grp_ppc(
  pp_site %>% filter(country == "NP"), "Sites New Zealand", "")

g1 <- (p1 | p2) / p3
g2 <- p4 / p5 / (p6 | p7)

pth1 <- file.path(
  "outputs", "figures", "outcomes", "secondary",
  "7-4-primary-model-acs-itt-ppc1.pdf")
pth2 <- file.path(
  "outputs", "figures", "outcomes", "secondary",
  "7-4-primary-model-acs-itt-ppc2.pdf")

ggsave(pth1, g1, width = 6, height = 5, device = cairo_pdf)
ggsave(pth2, g2, width = 6, height = 6, device = cairo_pdf)

g1

Code
g2

Sensitivity: Reduced Model (No site term)

The model is re-fit without site terms, but retaining the epoch terms.

Fit reduced model without site term
fit_no_site_model <- function(dat) {
  mdat <- make_primary_model_data(dat, c("inelgc3", "agegte60", "ctry"), c(10, 2.5, 1, 1))
  snk <- capture.output(
    mfit <- ordmod_epoch$sample(
      data = mdat,
      seed = 813508,
      chains = 8,
      parallel_chains = 8,
      iter_warmup = 1000,
      iter_sampling = 2500,
      refresh = 0,
      adapt_delta = 0.98
  ))
  mdrws <- as_draws_rvars(mfit$draws())
  
  # Add names
  epoch_map <- dat %>% dplyr::count(epoch, epoch_lab)
  names(mdrws$beta) <- colnames(mdat$X)
  names(mdrws$gamma_epoch) <- epoch_map$epoch_lab
  
  # Convert to treatment log odds ratios
  mdrws$Acon <- attr(mdat$X, "contrasts")$randA %**% mdrws$beta[grepl("randA[0-9]", names(mdrws$beta))]
  mdrws$Ccon <- attr(mdat$X, "contrasts")$randC %**% mdrws$beta[grepl("randC[0-9]", names(mdrws$beta))]
  mdrws$trtA <- mdrws$Acon[-1] - mdrws$Acon[1]
  mdrws$trtC <- mdrws$Ccon[-1] - mdrws$Ccon[1]
  
  mdrws$OR <- c(
    setNames(exp(mdrws$trtC), c("Intermediate", "Low with aspirin", "Therapeutic")),
    "Ineligible aspirin" = exp(mdrws$beta[which(grepl("inelg", colnames(mdat$X)))]),
    "Age \u2265 60" = exp(mdrws$beta[which(grepl("age", colnames(mdat$X)))]),
    setNames(exp(mdrws$beta[which(grepl("ctry", colnames(mdat$X)))]), c("AU/NZ", "NP"))
  )
  
  return(list(dat = mdat, fit = mfit, drws = mdrws))
}
res <- fit_no_site_model(acs_itt_nona_dat)
Odds ratio summary
odds_ratio_summary_table_rev(res$drws$OR)
Epoch and site terms
p <- plot_epoch_terms(exp(res$drws$gamma_epoch))
p

Sensitivity: Treatment coding

The model is re-fit with priors specified on the treatment contrasts rather than the orthonormal contrasts.

Fit primary model using treatment contrasts
res <- fit_primary_model(acs_itt_nona_dat, ctr = contr.treatment, save = FALSE)
Odds ratio table
odds_ratio_summary_table_rev(res$drws$OR)
Parameter Median 95% CrI Mean (SD) Pr(OR > 1)
Intermediate 1.15 (0.95, 1.40) 1.15 (0.12) 0.92
Low with aspirin 1.06 (0.81, 1.41) 1.07 (0.15) 0.67
Therapeutic 0.76 (0.42, 1.37) 0.79 (0.25) 0.18
Ineligible aspirin 1.08 (0.56, 2.08) 1.14 (0.39) 0.59
Age ≥ 60 0.58 (0.47, 0.72) 0.59 (0.06) 0.00
AU/NZ 0.76 (0.33, 1.76) 0.83 (0.38) 0.25
NP 0.68 (0.22, 2.55) 0.86 (0.70) 0.26
Posterior summaries for model parameters.
Epoch and site terms
p <- plot_epoch_site_terms(
  exp(res$drws$gamma_epoch),
  exp(res$drws$gamma_site),
  factor(res$dat$region_by_site, 
         labels = c("India", "Australia\nNew Zealand", "Nepal"))
)
p

Sensitivity: Concurrent controls

Three separate analyses are conducted based on concurrent randomisations as was done for the primary outcome. For these models, the epoch term has been removed.

Intermediate dose

In the following analyses, the set is restricted to participants who were randomised to either low-dose or intermediate-dose.

  • Set: ACS-ITT-intermediate
  • Covariates: anticoagulation intervention, age group, region of enrolment.
Descriptive summary table
save_tex_table(
  make_summary_table(acs_itt_concurc2_dat, "latex"), 
  file.path("outcomes", "secondary", "7-4-anticoagulation-concurrent-intermediate-summary"))
make_summary_table(acs_itt_concurc2_dat)
Anticoagulation intervention Patients Known Deaths DAFH, Median (Q1, Q3)
Low-dose 610 596 19 (3%) 23 (21, 24)
Intermediate-dose 613 603 15 (2%) 23 (21, 24)
Overall 1223 1199 34 (3%) 23 (21, 24)
Fit model
fit_acs_itt_intermediate <- function(fit = TRUE) {
  ctr <- contr.orthonorm
  X <- model.matrix( 
    ~ CAssignment + agegte60 + ctry, 
    data = acs_itt_concurc2_nona_dat, 
    contrasts.arg = list(CAssignment = ctr))[, -1]
  y_raw <- ordered(acs_itt_concurc2_nona_dat$out_dafh)
  y_mod <- as.integer(y_raw)
  N <- dim(X)[1]
  K <- dim(X)[2]
  unique_y <- length(unique(y_mod))
  mdat <- list(
    N = N, K = K, J = unique_y, X = X, y = y_mod, y_raw = y_raw,
    p_par = 2 * rep(1 / unique_y, unique_y),
    beta_sd = c(1, 2.5, 1, 1)
  )
  if (fit) {
    snk <- capture.output(
    mfit <- ordmod0$sample(
      data = mdat,
      seed = 75136,
      chains = 8,
      parallel_chains = min(8, parallel::detectCores()),
      iter_warmup = 1000,
      iter_sampling = 2500,
      refresh = 0
    ))
    mfit$save_object(file.path("outputs", "models", "secondary", "7-4-dafh", "acs-itt-intermediate.rds"))
  } else {
    mft <- readRDS(file.path("outputs", "models", "secondary", "7-4-dafh", "acs-itt-intermediate.rds"))
  }

  mdrws <- as_draws_rvars(mfit$draws())
  names(mdrws$beta) <- colnames(X)
  mdrws$Ccon <- as.vector(ctr(2) %**% mdrws$beta[1])
  mdrws$Ctrt <- mdrws$Ccon[-1] - mdrws$Ccon[1]
  mdrws$compare <- c("Intermediate vs low" = exp(mdrws$Ctrt[1]))
  mdrws$OR <- c(
    "Intermediate" = exp(mdrws$Ctrt), 
    setNames(exp(mdrws$beta[2:4]), c("Age \u2265 60", "Australia/New Zealand", "Nepal"))) 
  return(list(dat = mdat, fit = mfit, drws = mdrws))
}

res <- fit_acs_itt_intermediate(TRUE)
mdat <- res$dat
mfit <- res$fit
mdrws <- res$drws
Sample distribution
pdat <- acs_itt_concurc2_nona_dat %>%
  dplyr::count(
    CAssignment = factor(
      CAssignment, 
      levels = c("C1", "C2"),
      labels = c("Low", "Intermediate")),
    dafh = ordered(out_dafh, levels = 0:27), 
    .drop = F) %>%
  group_by(CAssignment) %>%
  mutate(p = n / sum(n)) %>%
  mutate(cp = cumsum(p)) %>%
  ungroup()
p <- ggplot(pdat, aes(CAssignment, p)) +
  geom_col(aes(fill = fct_rev(dafh)), colour = NA) +
  labs(x = "Anticoagulation intervention", y = "Cumulative proportion") +
  scale_fill_viridis_d("Days alive and\nfree of hospital", option = "A", direction = -1) +
  guides(fill = guide_legend(reverse = TRUE, ncol = 3)) +
  theme(legend.key.size = unit(0.75, "lines"))
p

Sample distribution
pth <- file.path("outputs", "figures", "outcomes", "secondary")
ggsave(file.path(pth, "outcome-7-4-descriptive-concurrent-intermediate.pdf"), p, height = 2.5, width = 4)
Posterior contrast
p <- plot_or_densities(mdrws$compare) +
  # geom_segment(aes(x = 1, xend = Inf, y = 1 - 0.1, yend = 1 - 0.1),
  #              arrow = arrow(length = unit(0.5, "lines"))) +
  geom_text(aes(x = 1, y = 1 - 0.1, label = sprintf("Pr(OR > 1) = %.2f", Pr(RV > 1))), hjust = 0, nudge_x = 0.01,
            family = "Palatino") +
  scale_x_log10(breaks = c(0.8, 0.9, 1, 1.1, 1.2, 1.3, 1.4, 1.5, 1.6)) +
  labs(y = "Comparison")
p

Odds ratio summary table
odds_ratio_summary_table_rev(mdrws$OR)
Parameter Median 95% CrI Mean (SD) Pr(OR > 1)
Intermediate 1.17 (0.95, 1.42) 1.17 (0.12) 0.93
Age ≥ 60 0.51 (0.41, 0.64) 0.52 (0.06) 0.00
Australia/New Zealand 0.61 (0.43, 0.87) 0.62 (0.11) 0.00
Nepal 0.60 (0.40, 0.91) 0.61 (0.13) 0.01
Posterior summaries for model parameters.
Odds ratio summary table
save_tex_table(
  odds_ratio_summary_table_rev(mdrws$OR, "latex"),
  "outcomes/secondary/7-4-primary-model-acs-itt-concurrent-intermediate-summary-table")
Code
mfit$summary(c("beta"))
# A tibble: 4 × 10
  variable   mean median     sd    mad      q5    q95  rhat ess_bulk ess_tail
  <chr>     <dbl>  <dbl>  <dbl>  <dbl>   <dbl>  <dbl> <dbl>    <dbl>    <dbl>
1 beta[1]   0.108  0.108 0.0724 0.0732 -0.0113  0.227  1.00   23232.   15548.
2 beta[2]  -0.666 -0.667 0.116  0.117  -0.858  -0.479  1.00   15313.   14791.
3 beta[3]  -0.487 -0.488 0.177  0.175  -0.778  -0.193  1.00   21270.   15933.
4 beta[4]  -0.509 -0.511 0.210  0.211  -0.850  -0.158  1.00   21021.   15504.
Code
mfit$diagnostic_summary()
$num_divergent
[1] 0 0 0 0 0 0 0 0

$num_max_treedepth
[1] 0 0 0 0 0 0 0 0

$ebfmi
[1] 0.9410 0.9585 0.9558 0.9978 0.9629 1.0617 0.9636 1.0391
Code
mcmc_trace(mdrws["beta"])

Code
mcmc_trace(mdrws["alpha"])

Low-dose with aspirin

  • Set: ACS-ITT-aspirin
  • Covariates: anticoagulation intervention, age group, region of enrolment.
Descriptive summary table
save_tex_table(
  make_summary_table(acs_itt_concurc3_dat, "latex"), 
  file.path("outcomes", "secondary", "7-4-anticoagulation-concurrent-stdaspirin-summary"))
make_summary_table(acs_itt_concurc3_dat)
Anticoagulation intervention Patients Known Deaths DAFH, Median (Q1, Q3)
Low-dose 299 291 11 (4%) 22 (20, 24)
Intermediate-dose 298 293 12 (4%) 22 (21, 24)
Low-dose with aspirin 283 280 10 (4%) 22 (20, 24)
Overall 880 864 33 (4%) 22 (20, 24)
Fit model
ctr <- contr.orthonorm
X <- model.matrix( 
  ~ CAssignment + agegte60 + ctry, 
  data = acs_itt_concurc3_nona_dat, 
  contrasts.arg = list(CAssignment = ctr))[, -1]
y_raw <- ordered(acs_itt_concurc3_nona_dat$out_dafh)
y_mod <- as.integer(y_raw)
N <- dim(X)[1]
K <- dim(X)[2]
unique_y <- length(unique(y_mod))
mdat <- list(
  N = N, K = K, J = unique_y, X = X, y = y_mod, y_raw = y_raw,
  p_par = 2 * rep(1 / unique_y, unique_y),
  beta_sd = c(1, 1, 2.5, 1)
)

snk <- capture.output(
  mfit <- ordmod0$sample(
    data = mdat,
    seed = 135356,
    chains = 8,
    parallel_chains = 8,
    iter_warmup = 1000,
    iter_sampling = 2500,
    refresh = 0
))
mdrws <- as_draws_rvars(mfit$draws())
names(mdrws$beta) <- colnames(X)

mdrws$Ccon <- as.vector(ctr(3) %**% mdrws$beta[1:2])
mdrws$Ctrt <- mdrws$Ccon[-1] - mdrws$Ccon[1]

mdrws$compare <- c("Intermediate vs low" = exp(mdrws$Ctrt[1]),
                   "Low with aspirin vs low" = exp(mdrws$Ctrt[2]),
                   "Intermediate vs low with aspirin" = exp(mdrws$Ctrt[1] - mdrws$Ctrt[2]))

mdrws$OR <- c("Intermediate" = exp(mdrws$Ctrt[1]),
              "Low with aspirin" = exp(mdrws$Ctrt[2]),
              setNames(exp(mdrws$beta[3:4]), c("Age \u2265 60", "Australia/New Zealand")))
Sample distribution
pdat <- acs_itt_concurc3_nona_dat %>%
  dplyr::count(
    CAssignment = factor(
      CAssignment, 
      levels = c("C1", "C2", "C3"),
      labels = c("Low", "Intermediate", "Low\nwith aspirin")),
    dafh = ordered(out_dafh, levels = 0:27), 
    .drop = F) %>%
  group_by(CAssignment) %>%
  mutate(p = n / sum(n)) %>%
  mutate(cp = cumsum(p)) %>%
  ungroup()
p <- ggplot(pdat, aes(CAssignment, p)) +
  geom_col(aes(fill = fct_rev(dafh))) +
  labs(x = "Anticoagulation intervention", y = "Cumulative proportion") +
  scale_fill_viridis_d("Days alive and\nfree of hospital", option = "A", direction = -1) +
  guides(fill = guide_legend(reverse = TRUE, ncol = 3)) +
  theme(legend.key.size = unit(0.75, "lines"))
p

Sample distribution
pth <- file.path("outputs", "figures", "outcomes", "secondary")
ggsave(file.path(pth, "outcome-7-4-descriptive-concurrent-stdaspirin.pdf"), p, height = 2.5, width = 6)
Posterior contrast
plot_or_densities(mdrws$compare)

Odds ratio summary table
odds_ratio_summary_table_rev(mdrws$OR)
Parameter Median 95% CrI Mean (SD) Pr(OR > 1)
Intermediate 1.30 (0.98, 1.72) 1.32 (0.19) 0.97
Low with aspirin 1.12 (0.84, 1.49) 1.13 (0.17) 0.78
Age ≥ 60 0.57 (0.44, 0.74) 0.57 (0.08) 0.00
Australia/New Zealand 0.97 (0.60, 1.57) 1.00 (0.25) 0.45
Posterior summaries for model parameters.
Odds ratio summary table
save_tex_table(
  odds_ratio_summary_table_rev(mdrws$OR, "latex"),
  "outcomes/secondary/7-4-primary-model-acs-itt-concurrent-stdaspirin-summary-table")
Code
mfit$summary(c("beta"))
# A tibble: 4 × 10
  variable    mean  median    sd   mad     q5     q95  rhat ess_bulk ess_tail
  <chr>      <dbl>   <dbl> <dbl> <dbl>  <dbl>   <dbl> <dbl>    <dbl>    <dbl>
1 beta[1]  -0.106  -0.105  0.104 0.104 -0.278  0.0644  1.00   24930.   15085.
2 beta[2]  -0.154  -0.155  0.101 0.100 -0.321  0.0129  1.00   26625.   15578.
3 beta[3]  -0.569  -0.569  0.136 0.139 -0.792 -0.344   1.00   17703.   15810.
4 beta[4]  -0.0354 -0.0349 0.248 0.248 -0.439  0.372   1.00   23944.   15070.
Code
mfit$diagnostic_summary()
$num_divergent
[1] 0 0 0 0 0 0 0 0

$num_max_treedepth
[1] 0 0 0 0 0 0 0 0

$ebfmi
[1] 1.0212 1.0008 1.0335 0.9707 1.0682 1.0523 1.0225 1.0058
Code
mcmc_trace(mdrws["beta"])

Code
mcmc_trace(mdrws["alpha"])

Therapeutic dose

  • Set: ACS-ITT-therapeutic
  • Covariates: anticoagulation intervention, age group, region of enrolment.
Descriptive summary table
save_tex_table(
  make_summary_table(acs_itt_concurc4_dat, "latex"), 
  file.path("outcomes", "secondary", "7-4-anticoagulation-concurrent-therapeutic-summary"))
make_summary_table(acs_itt_concurc4_dat)
Anticoagulation intervention Patients Known Deaths DAFH, Median (Q1, Q3)
Low-dose 79 75 3 (4%) 22 (20, 24)
Intermediate-dose 65 63 1 (2%) 22 (20, 24)
Therapeutic-dose 50 50 6 (12%) 22 (19, 24)
Overall 194 188 10 (5%) 22 (19, 24)
Fit model
ctr <- contr.orthonorm
X <- model.matrix( 
  ~ CAssignment + agegte60 + ctry, 
  data = acs_itt_concurc4_nona_dat, 
  contrasts.arg = list(CAssignment = ctr))[, -1]
y_raw <- ordered(acs_itt_concurc4_nona_dat$out_dafh)
y_mod <- as.integer(y_raw)
N <- dim(X)[1]
K <- dim(X)[2]
unique_y <- length(unique(y_mod))
mdat <- list(
  N = N, K = K, J = unique_y, X = X, y = y_mod, y_raw = y_raw,
  p_par = 2 * rep(1 / unique_y, unique_y),
  beta_sd = c(1, 1, 2.5, 1, 1)
)

snk <- capture.output(
  mfit <- ordmod0$sample(
    data = mdat,
    seed = 49135,
    chains = 8,
    parallel_chains = 8,
    iter_warmup = 1000,
    iter_sampling = 2500,
    refresh = 0
))
mdrws <- as_draws_rvars(mfit$draws())
names(mdrws$beta) <- colnames(X)

mdrws$Ccon <- as.vector(ctr(3) %**% mdrws$beta[1:2])
mdrws$Ctrt <- mdrws$Ccon[-1] - mdrws$Ccon[1]

mdrws$compare <- c("Intermediate vs low" = exp(mdrws$Ctrt[1]),
                   "Therapeutic vs low" = exp(mdrws$Ctrt[2]),
                   "Intermediate vs therapeutic" = exp(mdrws$Ctrt[1] - mdrws$Ctrt[2]))

mdrws$OR <- c("Intermediate" = exp(mdrws$Ctrt[1]),
              "Therapeutic" = exp(mdrws$Ctrt[2]),
              setNames(exp(mdrws$beta[3:5]), c("Age \u2265 60", "India", "Australia/New Zealand")))
Sample distribution
pdat <- acs_itt_concurc4_nona_dat %>%
  dplyr::count(
    CAssignment = factor(
      CAssignment, 
      levels = c("C1", "C2", "C4"),
      labels = c("low", "Intermediate", "Therapeutic")),
    dafh = ordered(out_dafh, levels = 0:27), 
    .drop = F) %>%
  group_by(CAssignment) %>%
  mutate(p = n / sum(n)) %>%
  mutate(cp = cumsum(p)) %>%
  ungroup()
p <- ggplot(pdat, aes(CAssignment, p)) +
  geom_col(aes(fill = fct_rev(dafh))) +
  labs(x = "Anticoagulation intervention", y = "Cumulative proportion") +
  scale_fill_viridis_d("Days alive and\nfree of hospital", option = "A", direction = -1) +
  guides(fill = guide_legend(reverse = TRUE, ncol = 3)) +
  theme(legend.key.size = unit(0.75, "lines"))
pth <- file.path("outputs", "figures", "outcomes", "secondary")
ggsave(file.path(pth, "outcome-7-4-descriptive-concurrent-therapeutic.pdf"), p, height = 2.5, width = 6)
p

Posterior contrasts
plot_or_densities(mdrws$compare)

Odds ratio summary table
save_tex_table(
  odds_ratio_summary_table_rev(mdrws$OR, "latex"),
  "outcomes/secondary/7-4-primary-model-acs-itt-concurrent-therapeutic-summary-table")
odds_ratio_summary_table_rev(mdrws$OR)
Parameter Median 95% CrI Mean (SD) Pr(OR > 1)
Intermediate 1.17 (0.65, 2.11) 1.22 (0.38) 0.70
Therapeutic 0.85 (0.46, 1.58) 0.90 (0.29) 0.31
Age ≥ 60 0.44 (0.25, 0.75) 0.46 (0.13) 0.00
India 1.60 (0.74, 3.52) 1.73 (0.72) 0.88
Australia/New Zealand 1.30 (0.75, 2.25) 1.35 (0.39) 0.82
Posterior summaries for model parameters.
Code
mfit$summary(c("beta"))
# A tibble: 5 × 10
  variable      mean    median    sd   mad     q5    q95  rhat ess_bulk ess_tail
  <chr>        <dbl>     <dbl> <dbl> <dbl>  <dbl>  <dbl> <dbl>    <dbl>    <dbl>
1 beta[1]  -0.222    -0.222    0.234 0.234 -0.609  0.162  1.00   23418.   15560.
2 beta[2]   0.000513  0.000222 0.213 0.212 -0.352  0.349  1.00   19752.   14656.
3 beta[3]  -0.820    -0.818    0.276 0.276 -1.28  -0.374  1.00   14620.   15492.
4 beta[4]   0.471     0.467    0.398 0.398 -0.182  1.13   1.00   18809.   14749.
5 beta[5]   0.260     0.259    0.282 0.280 -0.202  0.722  1.00   16305.   15958.
Code
mfit$diagnostic_summary()
$num_divergent
[1] 0 0 0 0 0 0 0 0

$num_max_treedepth
[1] 0 0 0 0 0 0 0 0

$ebfmi
[1] 0.9602 0.9842 0.9251 0.9629 1.0004 0.9787 0.9783 0.9464
Code
mcmc_trace(mdrws["beta"])

Code
mcmc_trace(mdrws["alpha"])